Ground states and formal duality relations in the Gaussian core model 
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We study dimensional trends in ground states for soft-matter systems. Specifically, using a high-dimensional 
version of Parrinello-Rahman dynamics, we investigate the behavior of the Gaussian core model in up to eight 
dimensions. The results include unexpected geometric structures, with surprising anisotropy as well as formal 
duality relations. These duality relations suggest that the Gaussian core model possesses unexplored symmetries, 
and they have implications for a broad range of soft-core potentials. 
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I. INTRODUCTION 



Soft-matter systems are notoriously difficult to analyze the- 
oretically, and much of what we know about their phase di- 
agrams is based on numerical simulations. Even classical 
ground states can almost never be derived from first princi- 
ples (see Refs. QjJ— J3[] for some rare exceptions). In this article, 
we place phenomena such as crystallization and solid-solid 
phase transitions in a broader context by studying dimensional 
trends in the Gaussian core model [6], in which particles in- 
teract via a Gaussian pair potential. 

The Gaussian potential models the entropic effective inter- 
action between the centers of mass of polymers |7[], and it is 
one of the simplest and most elegant soft-core potentials. The 
behavior of the Gaussian core model in two and three dimen- 
sions is relatively well understood (see, for example, Refs. iH] 
and [9]), but in higher dimensions, it remains mysterious. Di- 
mensions above three are an excellent test case for the study 
of phenomena such as decorrelation jioll , and this fits into the 
long tradition in statistical mechanics of studying the effect of 
dimensionality in interacting systems, such as critical dimen- 
sions for mean-field behavior (see Section 16.7 in Ref. ifTTll 
for an overview). 

Furthermore, higher dimensions play a fundamental role in 
information theory. Sphere packing is the low-density limit- 
ing case of the Gaussian core model, and sphere packings are 
error-correcting codes for a continuous communication chan- 
nel. The dimension of the ambient space for the packing de- 
pends on the channel and coding method used, and it can be 
quite high in practice fl2ll (up to thousands of dimensions). 
Thus, coding theory is a powerful motivation for the study of 
the high-dimensional Gaussian core model. 

Our conclusions are based on molecular dynamics simu- 



lations 1 13]. Such simulations are frequently used and often 
highly informative, but the computational difficulties are im- 
mense for many-body systems. Thanks to the curse of di- 
mensionality IU4I1 . high-dimensional simulations typically re- 
quire exponentially many particles, which severely limits the 
range of dimensions in which simulations are possible. To 
address this problem, we use a high-dimensional version of 
Parrinello-Rahman dynamics lfl5l fl6ll . Instead of imposing 
periodic boundary conditions using a fixed background lattice, 
we dynamically update the lattice using the intrinsic geome- 
try of the space of lattices. By increasing the adaptivity of the 
simulation, we are able to minimize the number of particles 
and avoid unnecessary computational complexity. This lets 
us carry out higher-dimensional simulations than were previ- 
ously possible. 

In this article we carry out Parrinello-Rahman simulations 
of the Gaussian core model in dimensions two through eight. 
In addition to observing surprising geometrical phenomena 
such as anisotropy, we find formal duality relations between 
Gaussian core ground states at densities p and 1 / p . Although 
such duality is known between reciprocal Bravais lattices (see, 
for example, Ref. lfl7ll ). it rarely holds for other structures, and 
its occurrence here suggests a deeper, not yet understood sym- 
metry of the Gaussian core model itself. 

Our approach fits into a program pioneered by Gottwald et 
al. 1U8I1 and applied in Refs. [ 19j,|20|]. They use genetic algo- 
rithms to search the space of candidate structures into which 
a fluid can freeze. Our goals are similar, but we make use 
of more analytic tools. Specifically, we compute gradients in 
the space of lattices, which enables us to use more powerful 
optimization techniques such as gradient descent or conjugate 
gradient. 



II. FRAMEWORK 
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Consider a periodic configuration of particles in n- 
dimensional Euclidean space R". Such a configuration is 
specified by an underlying Bravais lattice A C M", together 
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with a collection of translation vectors Vi,... , v#. In crystal- 
lographic terms, it is a lattice with basis. The particles are 
located at the points x + Vj for x G A and 1 < i < N. 

Given a radial pair potential V, the average energy per par- 
ticle is 

SfEI I VQz+v t -vj\). 
ZJV i=l ;'=1 zeA 

z^Oif i = j 

This quantity is the potential energy of the system, and the 
configuration is a classical ground state if it minimizes poten- 
tial energy, even allowing A and N to vary but keeping the 
particle density fixed. [The density equals N / det(A), where 
det(A) is the absolute value of the determinant of a basis for 
A.] In other words, classical ground states correspond to the 
canonical ensemble at zero temperature. 

Simulations typically fix A and allow v\,...,Vn to vary. 
This amounts to using A to define periodic boundary condi- 
tions. The lattice A remains the same throughout this process, 
and it is often chosen to be proportional to a hypercubic lat- 
tice IP for computational simplicity. This imposes artificial 
structure on the system, and N must be chosen quite large to 
minimize the effects of this structure. For example, if one 
wants the lattice spacing in A to be an order of magnitude 
larger than the typical spacing between particles, then N must 
grow roughly like 10". For large n, this is clearly infeasible, 
and even for n = 6, it requires careful use of all available com- 
putational improvements. Many simulations, such as those in 
Ref. yJJ , are therefore limited to roughly six dimensions. 

By allowing A to vary, one might hope to use a much 
smaller value of N. In the most extreme case, one could 
take N = I and study all Bravais lattice configurations. The 
naive dynamics then fail completely, because the forces in a 
Bravais lattice balance perfectly. Nevertheless, the potential 
energy varies dramatically between different Bravais lattices, 
with corresponding dynamics on the space of lattices. The 
Parrinello-Rahman method uses these dynamics. It can there- 
fore simultaneously update the underlying Bravais lattice A 
and the particle locations Vi , . . . , Vjy. 

Before describing the simulation results, we will give a 
derivation of the n-dimensional version of Parrinello-Rahman 
dynamics. It is equivalent to the original formulation in 
Refs. lfl5l[l6ll . except of course for the change in dimension. 
Instead of deriving it from a postulated Lagrangian, we show 
how it follows naturally from the intrinsic geometry of the 
space of lattices. Presenting the derivation gives us an op- 
portunity to describe some of the computational issues that 
become important in higher dimensions, such as the use of 
lattice basis reduction algorithms. 

HI. GEOMETRY OF THE SPACE OF LATTICES 

We will represent Bravais lattices by positive-definite, sym- 
metric matrices. Specifically, there is a linear transformation 
with matrix T such that A = TIP, If v = Tw, then the squared 
vector length v'v (the exponent t denotes transpose and we use 
column vectors) is w'T'Tw. Set G = T'T. This Gram matrix 



represents the metric in coordinates in which the underlying 
lattice is TP. Instead of fixing the metric and deforming the 
lattice, we will fix the lattice and deform the metric. This ap- 
proach is used to define the intrinsic geometry on the space of 
lattices (see, for example, Ref. [22]), and it makes the formu- 
las quite a bit simpler. 

To simplify the notation, define the function / of squared 
distance by /(*) = V(y/s)/2. Furthermore, write the vectors 
v\, . . . , vjv in the new coordinates as v, = Jw,-. Now the poten- 
tial energy of the system is 

w 

where we sum over all vectors w of the form w/ — uj +x with 
1 < iJ<N,xeZ n ,mdx^Oifi = j. 

The gradient of this sum as a function of G equals the ma- 
trix 

VC/(G) = ^Y i f(w'Gw)ww t . 

To see why, note that if we vary the i,j component G, ; 
of G while leaving all other entries fixed [and write w — 
{w\ , . . . , w„)], we find that 

f(w'Gw) = f \W Gw)wiWj. 

dGij 

When we update the configuration, we must fix det(G), 
so that the density of the system does not change [note that 
det(G) = det(A) 2 ]. However, the gradient does not respect 
this constraint. Instead, we must use the modified gradient 
VU(G) conditioned on fixing the determinant, which is com- 
puted as follows. If we define the standard inner product (•, ■) 
on the space of symmetric matrices by (A,B) = tr(AB), then 
to preserve det(G) we must remove the component of VU (G) 
in the direction of G 1 , because 

det(G + eVf/(G)) = det(G)(l +e(G~\ VC/(G)) + 0(e 2 )). 
We define the modified gradient VU (G) by 

so that (G- 1 ,Vf/(G)) =0. 

We could avoid this last complication by replacing the 
canonical ensemble with the grand canonical ensemble and 
controlling the particle density via the chemical potential. 
However, the modified gradient is not difficult to use, and our 
approach is convenient if one wishes to target a specific den- 
sity. 

IV. PARRINELLO-RAHMAN DYNAMICS 

The lattice gradient has two components, one for chang- 
ing G (computed above as the modified gradient) and one for 
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changing u\ , . . . , u^. For the latter, if we write w — ui — u> +x, 
then we need the gradient of f(w'Gw) as a function of iij and 
Ui, which is easily computed as follows. With w — u\ — ui+x, 
the gradient of f(w'Gw) as a function of m,- is 2f'(w'Gw)Gw. 
Thus, the w, component of the gradient of potential energy 
is the sum of (2/N)f / (w'Gw)Gw over all w of the form 
Ui — iij +x for some j plus the sum of — (2/N)f'(w'Gw)Gw 
over all w of the form uj — iij+x for some j. The full lattice 
gradient is made up of both the G and the u\,...,un compo- 
nents. 

Parrinello-Rahman dynamics consists of using the lattice 
gradient to define forces on the configuration. Thus, G 
changes as well as u\, . . . , w#. The simplest application is gra- 
dient descent, where we seek a local minimum for energy by 
following the negative gradient (or, better yet, using the con- 
jugate gradient algorithm). Of course, one could also use the 
forces in the usual way to define accelerations rather than sim- 
ply velocities, but we will focus on gradient descent here be- 
cause of our interest in ground states. 

The infinite sums in the algorithm must be truncated in 
practice, by summing only over the w such that WGw is 
at most some bound (chosen based on the decay rate of f). 
Writing w — Ui — uj + x, we must enumerate all 16Z" with 
this property. To do so, we use the Fincke-Pohst algorithm 
l23tl . which is far more efficient than brute-force searches. 
By contrast, simply summing over all the vectors in a large 
box becomes exponentially inefficient in high dimensions. We 
also periodically apply the L 3 lattice basis reduction algorithm 
lf24ll . which changes basis so as to keep the entries of G small, 
and we renormalize G to maintain a constant determinant (so 
that small numerical errors do not accumulate). 

As the simulation progresses, the metric changes and the 
connection with the original coordinates is lost. Nevertheless, 
using the Cholesky decomposition [25], we can recover the 
Bravais lattice TU' with basis Tu\,..., Tun from the matrix 
G and the vectors u\, . . . ,un by finding T such that G = T'T. 

Because of the need to use tools from lattice basis reduc- 
tion theory and linear algebra, the Parrinello-Rahman method 
is not as simple to implement in high dimensions as more 
straightforward methods are. However, in compensation it 
adapts itself to the structure of the system being considered 
and can therefore provide improved results. Furthermore, it is 
compatible with other standard computational methods such 
as Ewald summation or fast multipole methods II 1311 . 



V. RESULTS FOR THE GAUSSIAN CORE MODEL 

For the Gaussian core model, we take V(r) = exp(— nr 2 ) 
and hence f(s) = exp(— ns)/2. The choice of the constant % 
amounts to fixing the length scale, and it is chosen to make V 
self-dual under the Fourier transform. Let p denote the parti- 
cle density. 

The ground states of this model have been thoroughly ex- 
amined in up to three dimensions, although except in R 1 , no 
proof is known [5J. In R 2 , at all densities the ground state is 
the triangular lattice A%. In M?, the face-centered cubic lat- 
tice D3 is optimal at low densities, and the reciprocal body- 



n Energy 


State 


n Energy 


State 


1 0.04321740.. 


Z 


5 0.17434205... 


pc 2 


2 0.07979763 . . 


A 2 


6 0.19437337... 


5* 6 (1.0525...) 


3 0.11576766.. 


pc. 


7 0.21222702... 




4 0.14288224.. 


D A 


8 0.22788144... 


Eh 



TABLE I: Lowest known energies in dimension n when p = 1. The 
third column specifies the putative ground state, with "pcj" and 
"pc 2 " standing for phase coexistence between D3 and and be- 
tween D+ ( 1 .99750 . . . ) and D+ (0.50062 . . . ) , respectively. 

centered cubic lattice D\ is optimal at high densities lErjll . 
The cross-over point is at p = 1, but in fact the Maxwell 
double-tangent construction (i.e., the convexity of potential 
energy as a function of 1 /p) leads to phase coexistence for 
0.99899854... < p < 1.00100312.... This appears to give a 
complete description of the phase transition from D3 to . 

Little is known in higher dimensions, despite the connec- 
tions with coding and information theory. Cohn and Kumar 
JH] conjectured that the E$ and Leech lattices are universally 
optimal when n = 8 or 24, respectively. (In other words, they 
are ground states for the Gaussian core model at all densities. 
As shown in Ref. @], this implies optimality for many other 
potentials, such as all inverse power laws.) Torquato and Still- 
inger IU7I1 conjectured that in at most eight dimensions, certain 
Bravais lattices are always optimal at sufficiently high or low 
densities, but their conjecture was disproved in five and seven 
dimensions 112711 . Despite extensive exploration Il28ll . the true 
ground states have remained a mystery. Because of the dif- 
ficulty of simulation, previous studies have made use only of 
structures already known for other reasons. Comparing such 
structures in the Gaussian core model is of course of value, 
and it can sometimes lead to surprising results, but it provides 
little evidence as to the true ground states. 

We have run numerous Parrinello-Rahman simulations with 
2 < n < 8, 1 < < 24 (and occasionally larger), and various 
densities, with the following results (see also Table[l]». 

In two and three dimensions, we observe the previously 
known ground states. In four dimensions, we find the D4 
lattice at all densities, in accordance with the conjecture in 
Ref. lfl7ll . Thus, it appears probable that, like E$ and the Leech 
lattice, the D4 lattice is universally optimal. 

In five dimensions, we find different structures. The A 2 lat- 
tice was used in Ref. 1I27I1 to improve on Bravais lattices at 
low density; it consists of parallel translates of D4, repeating 
with period 4. Parrinello-Rahman simulations rapidly iden- 
tify and improve on this structure. It can be deformed by 
compressing the spacing between the parallel copies by some 
factor. A local minimum for energy is achieved for care- 
fully optimized values of the compression factor, which are 
between 0.998749 . . . and 1 when p < 1 and between 0.25 
and 0.250312 . . . when p > 1. Our simulations suggest that 
these are the true ground states, with the exception of phase 
coexistence for 0.99836946 . . . < p < 1 .00163526 . . .. 

These structures fit into the following general family. Let 

D„ = {(xi,. . . ,x n ) € Z" : x\ H Yx n is even} 
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FIG. 1: A plot of p- y l 2 {2U p + 1) as a function of lnp, where Up 
is the minimal energy attained by the structures D^(a) at density p. 
The reflection symmetry follows from formal duality. 

denote the checkerboard lattice in M", and let D+ be the union 
of D n with its translation by (1 /2, 1/2, . . ., 1 /2). Let D+(a) 
be D+ with the last coordinate scaled by a factor of a. That 
is, 

D+(oc) = {(xi,...,x„_i,ax„) : (*!,..., * n ) Gfl„ + }. 

Then (a) is the deformation of Aj with compression factor 
a/2. This is not obvious, but it can be checked by a straight- 
forward computation, and in fact it gives a substantially sim- 
pler construction of A5 than was previously known [namely, 
asD+(2)]. 

One noteworthy aspect of these configurations is their 
anisotropy. As the density increases, they experience greater 
compression along a distinguished axis than orthogonally to 
it. However, an even more surprising phenomenon is that 
there are formal duality relationships between these struc- 
tures. Formal duality is a generalization of the relationship 
between a Bravais lattice and its reciprocal lattice (see Sec.lVIl 
for more details). Formal duality relates the energies at densi- 
ties p and 1 / p : if E p is the Gaussian core energy at density p 
for a given structure and E p is that for its formal dual, then 

=p. 

2£ 1/p + l 

More generally, formal duality relates the energy of one struc- 
ture under a given pair potential to that of the formally dual 
structure under the Fourier transform of the potential. 

A non-Bravais lattice typically has no formal dual. Thus, it 
is remarkable that D+(a) is formally dual to Z)+(l / a). (See 
Sec. [VI] for a proof.) Formal duality implies that the high- 
density ground states in R 5 tend to D,j~(l/2), because the low- 
density ones tend to D<~(2) = Aj. FigureQ]illustrates the for- 
mal dualities within the D^(a) family of structures. 

In six dimensions, the lowest-energy states previously 
known were the Bravais lattice £g for low density and its 



reciprocal lattice E£ for high density, with a narrow region 
of phase coexistence in between. The lattices remain the 
lowest-energy states known at extreme densities, but in be- 
tween them, our simulations have identified other candidate 
ground states. They are all deformations of the orthogonal di- 
rect sum D3 D3 together with its translates by three vectors, 
namely, 

/l 1 1 1 1 1\ 

V2'2'2'2'2'2j ' 

(l,l,l,-i-i4)' and 

(These vectors are made up of holes in D3.) Define 3^^{a) to 
be the structure obtained by scaling the first three coordinates 
by a factor of a and the last three by 1 /a, so that volume is 
preserved. For certain values of a, these structures improve 
on E 6 and E* 6 for 0.25384516 . . . < p < 3.93940925 . . ., with 
phase coexistence for 3.93255017... < p < 3.94624440... 
and for the reciprocal range of densities. Unlike D^(a), 
which can be viewed as a modification of or Aj, the struc- 
tures ^(0:) differ more substantially from previously ana- 
lyzed structures. As in five dimensions, however, there are 
formal duality relations. Specifically, 3P(,{a) is formally dual 
to ^(l/a) for each a [of course ^(l/a) is isometric to 
^e(a)]. 

In seven dimensions, we find the D^(a) family of 
structures at all densities. For 0.04660088... < p < 
21 .45881937 . . ., the ground state seems to be itself (i.e., 
a = 1). For lower densities, we have a > 1 and for higher 
densities we have a < 1. Unlike the case of R 5 , there is no 
phase coexistence, because the optimal value of a changes 
continuously as a function of density. The low-density limit 
of D^(a) is Dt (V2), which is the same as the A^ structure 
studied in Ref . 12711 . By formal duality, the high-density limit 
isD+(l/V2). 

Finally, in eight dimensions our simulations provide further 
evidence that Eg (i.e., D%) is universally optimal. Parrinello- 
Rahman simulations are by no means limited to eight dimen- 
sions, and in fact, we have numerical results in as many as 
twelve dimensions. These results, together with a more ex- 
tensive analysis of the structures presented here, will appear 
elsewhere 02911 . 

Within the Z)+(a) family of structures for 1 < n < 8, the 
best energy at low density (i.e., the best sphere packing) is ob- 
tained when a = \J9 — n. For n < 4, the D+(\/9 — n) configu- 
ration is inferior to previously known sphere packings. How- 
ever, for 5 < n < 8 it achieves the highest sphere packing den- 
sity currently known. Note that D^(y/3) is the Ag packing, 
which has the same packing density as E& but is slightly infe- 
rior in the Gaussian core model at low densities. We have no 
conceptual explanation for why the six-dimensional behavior 
is subtly different from that in five, seven, or eight dimensions. 

It is also interesting to examine the a = 1 case. It seems 
that the and structures are not local optima for the 
Gaussian core model at any density. By contrast, Dj appears 
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to be the ground state over a large range of densities, and D% 
is almost certainly the ground state at all densities. It is pos- 
sible that is also universally optimal: so far we have not 
explored this case as thoroughly as those in lower dimensions, 
but we have not yet found any structure that beats D$ at any 
density. We will examine this issue elsewhere ll29ll . Univer- 
sal optimality cannot hold for D+ with n > 10, because these 
packings are not even optimal sphere packings. 



VI. FORMAL DUALITY 

Recall that Poisson summation relates the sum of a func- 
tion over a Bravais lattice to the sum of its Fourier transform 
over the reciprocal lattice. Specifically, given a sufficiently 
well-behaved function / : R" — > K (for example, a Schwartz 
function) and a Bravais lattice A C R", 



E/to 

xeA 



vol (I 



Here, vol(R"/A) denotes the volume of a fundamental do- 
main of A and we normalize the Fourier transform and recip- 
rocal lattice by 

760= / f(x)e 2 ^dx 



and 



A*={y£R": (x,y) £ Z for all x e A}. 

Formal duality is a generalization of Poisson summation to 
certain special structures that are not Bravais lattices. To state 
it correctly, it is important to view the sum 



I /to 



jeGA 

not simply as a sum over points in A, but rather as a sum 
over vectors between points in A. For a Bravais lattice these 
notions are exactly the same, but in more general cases they 
are not. (In fact, Poisson summation cannot be generalized 
from the first point of view BOH .) For example, consider a 
periodic configuration given by the union of N disjoint lattice 
translates A + Vj, . . . , A + v^. Then the analog of summing 
over the lattice is 

1 N N 

^ELE/(*+vy-v*). 



j=lk=lxeA 



Call this sum the average pair sum of / over the configura- 
tion. It is the average over all points in the configuration of 
the sum of / over all vectors from it to other points. (Thus, 
it is independent of how the configuration is decomposed into 
translates of Bravais lattices.) 

Suppose & and £2 are particle arrangements, where 2? has 
particle density 8 and £2 has particle density 1/5. We say & 
and £} are formal duals if for every Schwartz function /, the 
average pair sum of / over 9* is 8 times that for / over £}. 



Poisson summation shows that this definition generalizes the 
case of reciprocal Bravais lattices. If is formally dual to 
itself, we call it formally self-dual, and if it is formally dual 
to an isometric copy of itself, we call it formally isodual. For 
example, the triangular lattice in the plane is not self-dual, 
since its reciprocal lattice is a rotated copy of itself, but it is 
isodual. 

For any periodic configuration, one can always write the 
average pair sum of / in terms of /, by using a generalized 
Poisson summation formula: for a Bravais lattice A and trans- 
lation vector v, 



I /(a 

xeA 



-V) 



1 



vol(R«/A) 



(In fact, the right side is the Fourier expansion of the left side 
as a function of v that is periodic modulo A.) The average pair 
sum of / over A+ vi, . . . ,A + Vn then becomes 



N 



vol(KVA) y S/ () ' ) 



N 



The factor of N / vol(M"/A) is the density of the configuration, 
so the question becomes whether the remaining sum is the 
average pair sum of / over some periodic structure. Except 
when N = 1, it usually is not: for example, the coefficient 
of /(y) is generally irrational. Formal duality only arises in 
exceptional cases. It is not obvious when a configuration has 
a formal dual or, if it does have one, what the formal dual is. 

Proposition 1. The D+ structure is formally self-dual when n 
is odd or a multiple of four. When n is even but not a multiple 
of four, DJ is formally isodual. 

When n is even, D+ is a Bravais lattice, whose reciprocal 
lattice is D+ if n is a multiple of four and D+ (—1) otherwise. 
When n is odd, D+ is not a Bravais lattice and the duality is 
more subtle. 

Proof. Let v = (1/2, 1 /2, . . . , 1 /2), so Z)+ is the union of D„ 
and D„ + v. Then the average pair sum of / over D+ is 



7 E 760 i+e 2 ^ 2 = E 760- 



-cos(27r(v,y)) 



The reciprocal lattice D* consists of four translates of D„, by 
the vectors 0, v, (0,0, ... ,0, 1), and (1 /2, 1/2, . . . , 1 /2, -1/2). 
In each of these four cases, the inner product (v,y) is easily 
understood. In the first case, it is an integer, in the second it is 
an integer plus n /4, in the third it is an integer plus 1 /2, and 
in the fourth it is an integer plus (n — 2)/4. 

When n is odd, this yields weights of 1, 1/2, 0, and 1/2 
multiplying f(y) in the four cases. Because n is odd, the 
vector (1,1,. .. ,1,0) is in D n and hence translating D„ by 
(1/2, 1/2, 1/2, — 1/2) is equivalent to translating it by — v. 
Thus, the average pair sum for / is simply 

\ E (27to+7(*+v)+7(*-v)), 
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which is the same as the average pair sum for / over D+. It 
follows that D+ is formally self-dual when n is odd. 

An analogous computation works in the even case, or it can 
be verified more simply using the fact that is then a Bra- 
vais lattice. □ 

Lemma 2. If 8? and £1 are formally dual structures in M", 
and T : W — > W is an invertible linear transformation, then 
T and (T')~ l £} are formally dual. 

Here T' denotes the adjoint operator with respect to the in- 
ner product (i.e., the transposed matrix). 

Proof. To compute the average pair sum for / on T simply 
compute it for the composition fo T on The Fourier trans- 
form of fo T is (detr) -1 (/o (r') _1 ). Now applying formal 
duality for & and £2 completes the proof. □ 

It follows immediately that D+(a) is formally dual to 
(1 /a) when n is odd or a multiple of four: apply Lemma[2] 
with = J2 = D+ and with T being the map that multiplies 
the last coordinate by a. 

The formal isoduality of 3^(,{a) is proved similarly. For 
a = 1, formal self-duality follows from a calculation much 
like the proof for Proposition [T] Then Lemma [2] implies that 
&(,(a) is formally dual to SP^iX /a), but of course the two 
configurations are isometric. 

In the literature, formal duality is usually understood to re- 
fer only to radial functions / (see, for example, p. 185 of 
Ref. 13111 ). That is a weaker condition, which depends only 
on the radial pair correlation functions of the structures. We 
have defined a stronger version of formal duality in this paper, 
without that restriction, because the stronger version in fact 
holds for the structures we find in our simulations. (Further- 
more, it behaves better. For example, the proof of Lemma [2] 
breaks down in the radial case, because / and fo T will gener- 
ally not both be radial.) However, for the discussion in Sec.lVl 
the pair potential is isotropic, so only the radial version of for- 
mal duality is needed. Note also that radial symmetry erases 
the distinction between formal self-duality and formal isodu- 
ality. 

VII. CONCLUSIONS AND DISCUSSION 

We have used Parrinello-Rahman dynamics to identify 
ground states in dimensions that were previously beyond the 
reach of simulation. This approach is effective because it 
adapts to whichever underlying Bravais lattice is most favor- 
able. That means it probably offers little advantage in detect- 
ing disordered states or even phase coexistence, but it is appro- 
priate whenever one anticipates a high degree of symmetry. 

The formal duality relations are the most noteworthy con- 
sequence of our simulations. Such relations occur only rarely 
for structures other than Bravais lattices, and it is far from ob- 
vious why they arise here. A remarkable possibility is that 
all periodic ground states of the Gaussian core model, in any 
dimension, may occur in formally dual pairs. If true, this hy- 
pothesis deserves a more conceptual explanation than a case- 
by-case calculation, and it suggests that the model possesses 



deeper symmetries than are currently understood. Even if it 
is false, there must be a reason why formal duality arises so 
frequently in low dimensions. 

The families of structures studied in this article have much 
broader applicability than just to the Gaussian core model. We 
believe that they will minimize many other repulsive potential 
functions, such as inverse power laws, although we have done 
relatively little experimentation in this direction. 

One reason for our focus on the Gaussian core model is that 
it is the natural setting for studying universal optimality JH. 
The known universal optima include some of the most sym- 
metrical and beautiful geometrical configurations, with con- 
nections to many other topics such as sporadic finite simple 
groups and exceptional Lie algebras. Relatively few universal 
optima are known and any new examples are of interest. In 
this article, we have described simulation evidence that D4 is 
likely universally optimal and that Dg may be. Both cases are 
surprising: the analog in spherical geometry of the universal 
optimality of D4 turned out to be false [32], and there have 
not even been any previous hints that Dg might be universally 
optimal. 

Our results also offer insight into the complexity of ground 
states. One measure of the complexity of a lattice is the num- 
ber of Bravais lattice translates required to generate it (in crys- 
tallographic terms, the minimal size of a particle basis). Bra- 
vais lattices have complexity 1, while disordered structures 
can be considered to have infinite complexity. The hexago- 
nal close-packing has complexity 2, while the structures intro- 
duced here have complexity 2 (in five and seven dimensions) 
and 4 (in six dimensions). 

Do the complexities of ground states grow with dimen- 
sion? The Torquato-Stillinger decorrelation conjecture ifloll 
suggests that they do grow and eventually become infinite. If 
so, how quickly do they grow? There is a striking example 
in ten dimensions (the Best packing 113 ill , which is the dens- 
est sphere packing known in R 10 , and which has complexity 
40), but other low-dimensional ground states for repulsive po- 
tentials that have been reported in the literature typically have 
much smaller complexity, with the exception of phase coexis- 
tence. 

In up to eight dimensions, our results for the Gaussian core 
model suggest that the ground states may indeed have low 
complexity for most densities. The structures we have identi- 
fied seem difficult to improve, even if we allow the algorithm 
the freedom of substantially higher complexity, and we sus- 
pect that they are the true ground states. Of course, we cannot 
rule out the possibility that extraordinarily high-complexity 
states offer tiny improvements, but we consider it unlikely. 

We conclude with a computational challenge regarding 
simulation in high dimensions. It is undoubtedly impossible 
to carry out effective simulations in extremely high dimen- 
sions, but where is the threshold for feasibility? For example, 
is it possible in 24 dimensions? Many remarkable phenomena 
in mathematics and information theory (such as the Leech lat- 
tice 113 111 ) occur there, and reliable simulation results would be 
very interesting. 
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